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Abstract 

A numerical investigation of grain-boundary grooving by means of a Level Set 
method is carried out. An idealized polygranular interconnect which consists of grains 
separated by parallel grain boundaries aligned normal to the average orientation of the 
surface is considered. The surface diffusion is the only physical mechanism assumed. 
The surface diffusion is driven by surface curvature gradients, and a fixed surface slope 
and zero atomic flux are assumed at the groove root. The corresponding mathematical 
system is an initial boundary value problem for a two-dimensional Hamilton-Jacobi 
type equation. The results obtained are in good agreement with both Mullins' ana- 
lytical "small slope" solution of the linearized problem |14j (for the case of an isolated 
grain boundary) and with solution for the periodic array of grain boundaries (Hackney 

0). 



1 Introduction 

This paper presents the results of our work on numerical modeling and simulation of 
grain-boundary (GB) grooving by surface diffusion. Our ultimate goal is to develop 
and test a fast numerical approach for the simulation of formation and propagation of 
groove-like defects in thin film interconnects used in microelectronics (ME). 



In modern ME industry, the quality and reliability of ME integrated circuits have 
become no less important than their performance. Some of the most vulnerable el- 
ements of ME circuits, susceptible to several types of mechanical failures, are the 
interconnects. These are metallic conductors which connect the active elements. 

The defects (due to the small cross-section, high current density, mechanical stresses 
and presence of GBs acting as fast diffusion pathways) lead to the loss of electrical and 
mechanical integrity, i.e. to line opens or shorts. Thus, such defects are one of the 
main reliability concerns in advanced integrated circuits. 



1.1 Mechanisms of Mechanical Failure in Interconnect Lines 

In this section we describe some basic failure mechanisms in interconnects and outline 
an appropriate physical model. 

Many properties of polycrystalline materials are affected by the intersection of GBs 
with external surfaces, especially in the presence of applied or internal fields. Com- 
mon examples are growth of GB grooves and cavities [[H], o], stress voiding [27] and 
electromigration [g, |l^, 17, ^TJ. 



In the absence of an external potential field, the GB atomic flux Iqb = and the 
corresponding groove profile evolves via surface diffusion under well-known conditions 



of scale and temperature (the so-called Muffins' problem [14|). Mass transport by 
surface diffusion is driven by the surface Laplacian of curvature. Essentially, for convex 
surfaces, matter flows from high-curvature regions, while for concave surfaces the flow is 
from low curvature regions. In order to solve surface-diffusion problems, four different 
approaches have been taken. We refer the interested reader to the article by Zhang 



and Schneibel [28], where these approaches are discussed and to the references therein. 

The physical origins of a GB flux may be gradients of the normal stress at grain 
boundaries Q and/or electromigration forces 0. GB grooving with a GB flux in 
real thin film interconnects is a complex problem. It requires sophisticated numerical 
modeling technique which can manage with such issues as aperiodic arrays of GBs, 
anisotropy of the surface tension, GB migration, formation of slits with a local steady- 
state shape in the near-tip region and bridging across the slits near their intersections 



with the surface left behind [17|. Level Set Method seems to be a good candidate for 
addressing the problems, however it has never been used yet to this aim. As the first 
step in application of LS Method to the problem of grooving with EM flux, we test 
in this paper LS Method over two simple -and already solved- grooving problems and 
compare the LS Method' results with those obtained previously in Jl4[ , [|. First is 
classical Mullins' problem (GB grooving controlled by surface diffusion in an infinite 
bicrystal with a stationary GB). Second is GB grooving by surface diffusion in the 
periodic GB array of stationary GBs. The electromigration flux will be taken into 
account in the next publication. 

Below we give more details related to the physical model. 

• Driving Forces and Diffusion Mobilities 

In the absence of an electric current, the diffusion is driven by a variation in 
chemical potential, fi, which causes atoms to migrate from high potential to low 



potential regions. It may be shown that [14] 



ti(K) = K 7 n, (i.i) 
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where K is the surface curvature, 7 is the surface tension, and VL is the atomic 
volume. Gradients of chemical potential are therefore associated with gradients 
of curvature. 

In interconnects, GBs represent numerous fast diffusion pathways with high dif- 
fusion coefficient, D. As a matter of fact, the bulk diffusion can be neglected 
The diffusion flux along the GB, Iqb, is given by 

Igb = ^V^, (1.2) 

where 5 ~ 10 -8 cm is the GB thickness, k the Boltzmann constant and T the 
absolute temperature. 

Let r be the tangential direction to the surface profile in 2D. If n = (n x ,n y ) is 
the unit vector normal to the surface or GB, then the following relations hold: 

v dK dK dK 

T = (n y ,-n x ), — = VK • r = — n y - — n x = K T . (1.3) 

The surface diffusion flux along the groove walls is given by the formula 

J = -B K T , (1.4) 

where 

is known as Mullins' constant. Note that J is proportional to the first directional 
derivative of the curvature. 

Boundary Conditions 

The boundary conditions at the groove root are dictated by the local equilibrium 
between the surface tension, 7, and the GB tension, 7^. In the symmetric case 
of the GB (x = 0) normal to an original {y = const.) flat surface, the angle of 
inclination of the right branch of the surface at the groove root with respect to 
the x axis is 6q = sin~ 1 ('ygi ) /2j) (see Fig. ||). 

The rapid establishment of the equilibrium angle between the GB and the surface 
by atomic migration in the vicinity of the intersection develops some curvature 
gradient at the adjacent surface and thus induces a surface diffusion flux along 
the groove wall in the direction away from the groove root, opposite to the groove 
extension direction. 

Other boundary conditions depend on the particular problem (presence or absence 
of GB flux, etc.) 
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Figure 1: Sketch of a GB groove, w denotes the half-width of the groove, and d denotes the 
depth. 



2 Mathematical Model 

2.1 The Conventional Approaches 

An adequate mathematical model which captures the above physical phenomena in 
interconnects was developed first by Mullins and extended by him and others 
||To| , p| , |i~5f . It describes the evolution of the groove shape, y(x, t), and has the form of 
a transport equation 



yt 



-J x = -B{(l + ylr l l 2 [{l + yl)-^yJ } . 

I L J x ) x 



(2.1) 



J and B are given in (|1.4j ) and (p..q). 

For an isolated GB at x = 0, the groove continues to develop because the material 
continues to move from the curved shoulder of the groove to the flat surface. The 
classical description is provided by an analytic solution (on the x > side) of the 
linearized version of the equation fl2.1| ) (the "small slope approximation", SSA). The 
linearized equation has the form [|l4[ 

Vt = ~B y xxxx , (2.2) 

subject to the initial condition 

y(x,0) = const, (2-3) 

and the boundary conditions 

y x (0, t) = tan <9 = m < 1, (2.4) 
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J(0,t) = y xxx (0,t) = 0, 
y(x — » oo,t) = const with all derivatives. 



The first condition in ( |2.4D is the small slope approximation itself. The second one 
reflects the absence of a GB flux Igb- The solution describes a profile with a constant 
shape whose size is increasing all the time. 

Although this analytical approach describes some basic phenomena in interconnects, 
it is of limited use due to the restriction on the steepness of the slope. There are 
several numerical techniques which are widely used in modeling moving fronts, such as 
the marker/string (M/S) methods [22] or the volume- of -fluid (VOF) methods || jlBf . 



These methods deal directly with the evolution equation of type ( |2.1| ), and therefore 
are "explicit" methods. The M/S methods come from Lagrangian approach to front 
evolution problems. In the Lagrangian approach, the grid is attached to the moving 
front. A known drawback of the Lagrangian approach is that it is not well-suited for 
the computation of bifurcating fronts. Besides, stability and local singularity problems 
are more emphasized in these methods than in methods based on Eulerian approach, 
such as the VOF method. The Eulerian approach, in which the front moves through a 
grid which is fixed in space, does not have these drawbacks, but - as it is known - here 
the fronts are diffused. In addition, some intricate (subcell) bookkeeping is required to 
properly keep track of fronts. 

There are numerical approaches which are based on finite- element discretization of 
the computational region Q. However, they result in complicated algorithms which 
involve many computational steps such as computations of the following: displacement 
field of material points from a reference configuration, the stress field as a result of 
diffusion in the solid and geometry update of interfaces. Besides, the computational 
complexity grows since higher resolution is required as the shape of the interface be- 
comes more complicated. As a result, these methods are unable to handle very complex 
multidimensional boundary shapes. 



2.2 The Proposed Solution: Usage of the Level Set Method 

To "capture" the interface (rather than to track it), our method of choice is the "im- 
plicit" LS method. The method was introduced by Osher and Sethian and was further 
developed during the last several years (for an introduction to the LS methods and 
an exhaustive bibliography list see the monographs by Sethian p3| , pi]]). The method 
enables to capture drastic changes in the shape of the curves (interfaces) and even 
topology changes. 

The basic idea of the method consists of embedding the curve y(x, t) into a higher 
dimensional space. As a matter of fact, we consider the evolution of a two-dimensional 
field (j)(x,y,t) such that its zero level set, (p(x,y,t) = 0, coincides with the curve 
of interest, y(x,t), at any time moment t. The level set function <p(x,y,t) can be 
interpreted as a signed distance from the curve y(x,t), which moves in the direction 
normal to itself. 

The evolution of 4>(x, y, t) is described by an Hamilton-Jacobi type equation. A 
remarkable trait of the method is that the function (p(x, y, t) remains smooth, while 
the level surface <f> = may change topology, break, merge, and form sharp corners as 
4> evolves. Thus, it is possible to perform numerical simulation on a discrete grid in 
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the spatial domain, and substitute finite difference approximation for the spatial and 
temporal derivatives in time and space. 
The evolution equation has the form 

</> t + F|Vc/>| = 0, given (j>(x,t = 0). (2.5) 

The normal velocity, F, is considered to be a function of spatial derivatives of 4>(x, y, t). 
In many applications F is a function of the curvature, K, and its spatial derivatives. 
The curvature K may be computed via the level set function as follows: 

K = V-n, n = ^-=\ Yx . /0 , — — -pr I . (2.6) 

|V " l^ + ^) 1/2 (^) 1/2 

Here n is "normal vector" , and it coincides with the (previously introduced) unit normal 
to the surface, y(x,t), on the zero level set <f) = 0. Formulas (|2.6|) can be combined as 
follows 

„ „ V4> 4>Xx4> 2 y ~ 2(j) X (f>y(j) X y + <j)y y 4>l 

v '™r — — ' (2 ' 7> 



and the sign of K is chosen such that a sphere has a positive mean curvature equal to 
its radius. In the case of surface diffusion in 2D, 

F = -BK TT . (2.8) 

One drawback of the LS method stems from its computational expense. Its com- 
plexity seems to be as much as 0(n 2 ) operations per time step, which is more than 
any Lagrangian method which necessitates 0{n) operations per time step, where n is 
the number of grid points in the spatial direction. It is possible, however, to reduce 
the complexity of the LS method to O(n) using a local (another term is narrow band 
(tube)) approach |20[| . This is achieved by the construction of an adaptive mesh 
around the propagating interface. We distinguish between the "near field", which is a 
thin band of neighboring level sets around the propagating front, and the "far field" 
which contains the rest of the grid points. The evolution equation is solved only in the 
near field. The values of 4> at grid points in the far field are not updated at all. When 
the interface in motion reaches the edge of the narrow band, a new narrow band is built 
around the current interface position. Note that this could be done without interface 
reconstruction from the level set function (which requires some additional computa- 
tions). We just have to examine the shift in the sign of at grid points adjacent to 
the interface. The width of the narrow band is determined as a balance between the 
computation involved in the re-built and the calculations performed on far away points. 

In most of the applications of the LS method to date, the driving forces were 



proportional to the curvature (see [23, 24] for review and discussion). There are only few 
applications [l^] where the driving force is proportional to the second directional 
derivative of the curvature (in the 3D case, to the surface Laplacian of curvature which 
is constructed from the derivatives in each principal direction), which is the case for 
the normal velocity function ( |2.8| ). Therefore, the present materials science problem 
presents a rather new (from the mathematical point of view) application for the LS 
method. As pointed out in |J, "this is an intrinsically difficult problem for three 
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reasons. First, owing to the lack of a nice maximum principle, an embedded curve 
need not stay embedded, and this has significant implications in attempting to analyze 
motion which results in topological change. Second, the equations of motion contain 
a fourth derivative term, and hence are highly sensitive to errors. Third, this fourth 
derivative term leads to schemes with very small time steps." 

2.3 Computational Algorithm 

A typical computational domain is a rectangular box [0, 1%; 0, 12] of a material in 2D. 
The proposed computational algorithm consists of the following steps: 

BEGIN ALGORITHM 

1. Discretization - The entire computational region W is discretized using a uni- 

form grid Xi = iAx,yj = jAy,i = O...N,j = 0...M, where N and M are the 
number of grid points in x- and y- directions respectively. The functions are 
projected on this grid, so that 4>(x,y,t) = (f>i t j(t). 

2. Initialization - The initial interface, y(x,t = 0), is defined analytically, or as a 

set of points in W (the points lie on x = const grid lines, but not necessarily on 
y = const grid lines). In the latter case, we define a cubic spline £(x,t = 0) passing 
through these points in order to be able to perform further initializations. The 
function y(x,t = 0) needs not to be necessarily smooth (i.e., it may feature sharp 
corners, discontinuities, etc.), but, in our implementation it must be single- valued 
in order to make it possible to choose the sign of <f> (below). This is because we 
are only interested in the particular case of analyzing the motion of open curves 
which may be described by functions during the whole process of the evolution. 
We also define the near field and the far field. The width of the near field is 
usually 5 to 10 grid levels (points). 

In the region W, the level set function (p is initialized as an exact signed distance 
function to the initial interface (see Fig. ||), 

4>(x uyj ,t = 0) <0 if y j <y(x,t = 0) (2.9) 

<j){xi,y jy t = 0) = if yj = y(x,t = 0) 
<p(xi,yj,t = 0) > if yj>y(x,t = 0). 

@). Since (p(x,y,t = 0) is a signed distance function, then \V(f>(x,y,t = 0)| = 1. 

3. Compute normal vector components and curvature using formulas J2.6|) , Q2.7| ). 

The derivatives in ( |2.6[) , (|2.7| ) (as well as in other functions of x,y except the 
gradient term in the evolution equation itself, see step 6) are discretized using the 
standard second order accurate central difference approximations. Fourth-order 
accurate approximations were tested also but we did not observe any particular 
increase in the global accuracy of the calculations. In addition, in this case, 
the implementation of the boundary conditions with the level set function is 
problematic due to the use of a wide stencil. The time step also needs to be 
reduced in order to have stability. We find that the standard central difference 
scheme works well for us. 
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Figure 2: Computational domain. 



4. Compute first directional derivative of the curvature, K T , using the formula (1.3) 
and second directional derivative of the curvature, K TT , 



r . n[v7!/ , -K XX 4> 2 + 2K X y<j) X <t> y - Kyy4>l K {K X (j) X + Ky(f)y) 

K TT = V VA • r • r = To \ — ri5 1" 



-K XX (j) 2 +2K X y(j) X (j) y -Kyy(g 

V K[K T + K y (n x + n y ) - K x {n y - n x )\ . 



(2.10) 



We now have the normal velocity function fl2.8|) and the flux (|1.4| ). 

5. Choose time step. The CFL condition for the surface diffusion is 

Aii < min 4 (Aa;, Ay) /B. (2.11) 

The CFL condition for the Hamilton-Jacobi equation in updating the velocity is 

At 2 < min(Ax, Ay)/F max , (2.12) 

where F max is the largest magnitude of the normal velocity in the computational 
domain. The adaptive time step At is chosen as the smallest of the two. 

6. Compute backward and forward gradient functions; update (ft from the evolution 

equation using explicit time-stepping scheme. The solutions of equation fl2.5| ) are 
often only uniformly continuous with discontinuous derivatives, no matter how 
smooth the initial data is ]l8| , [HJ. Simple central differencing is not appropriate 
here to approximate the spatial derivatives in |V0|. Instead, we use Essentially 
Non-Oscillatory (ENO) type schemes for Hamilton-Jacobi equations as developed 
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in Ijlq |l9|, pB| . More precisely, we use second-order ENO scheme given explicitly 



in [29|. To update <p for one time step, the simplest method is to use Euler, i.e. 

^n+l = + AfL(0 n ), (2.13) 

where L((f>) is the spatial operator in (|2.5|). 

7. Update near field. Check the sign of <p at the grid points adjacent to the interface 
and compute the new locations of near field points. 

Go to step 3 

END ALGORITHM 

Remark 1: To achieve a uniformly high-order accuracy in time, we replace ( |2.13| ) with 
the second-order Total Variation Diminishing (TVD) Runge-Kutta type discretization 
fl?| p5| , which reads 



At r 



=( lu + 



2 

The necessary changes to the algorithm are obvious. The choice of such a low-order 
Runge-Kutta scheme is justified by the fact that the time step, dictated by stability 
requirements, is very small. 

Remark 2: It is highly desirable that the level sets behave nicely, in the sense that 
two different level sets do not cross, and in fact remain roughly evenly spaced in time. 
In terms of the level set function (f>, this corresponds to the fact that the gradient of 
cf) at any given point of a level set does not change dramatically over time. For the 
numerical method this translates into numerical stability. The best way to achieve this 
is to keep 4> close to the signed distance function (or even to keep it exactly equal to 
the signed distance function), thus keeping |V^| ~ (=)!• The operations performed 
on <j) that accomplish it are called "reinitialization". To summarize, reinitialization is 
the process of replacing 4>(x, y, t) by another function cp(x, y, t) that has the same zero 
contour as <p{x, y, t) but behaves better, and then taking this new function (j>(x, y, t) as 
the initial data to use until the next round of reinitialization. There are several ways 



to do this. The straightforward one (first proposed in [13] and recently used in [g]) 
is to interrupt the time-stepping, reconstruct the interface using some interpolation 
technique and directly compute a new signed distance function to the interface. This 
approach is very expensive and also may bring some undesirable side effects, such as 



oscillations in the curvature. Instead, we use the iteration procedure of [26|. The 
function <f> is reinitialized by solving the following Hamilton-Jacobi type equation to 
its steady state, which is the desired signed distance function: 

4> t = Sfa) (1 - |V^|) , (2.15) 

where S is a smoothed sign function 

5(0o) = Jtl , e = min(A 2 ;, Ay). (2.16) 

The same second-order ENO and TVD Runge-Kutta schemes used for the solution 
of the equation ( |2.5| ) are used for the iteration of ( |2.15| ), As a rule, three or four 
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Figure 3: The cosine curve, evolving under (2J5) with F = —0.1 K TT . A coarse 75x75 grid 



is used. 25000 time steps were made by the Runge-Kutta integrator (the shape is printed 
out every 500 steps), and we reinitialize in every step. 

iterations are sufficient to evolve 4> close enough to the desired signed distance function. 
An important practical question is how frequently the reinitializations are applied. In 
some applications of the level set method the reinitializations could be triggered after 
a fixed number of time steps. However, we achieved the best results by reinitializing 
every time step in the band of level sets that contains points from the near field. 
Remark 3: The evolving interface touches the vertical boundaries x = 0, x = l\ by 
its ends and therefore any boundary conditions imposed on vertical walls influence 
the evolution of the front. This is why, depending on the nature of the problem, we 
either choose periodic b.c. at vertical walls or just an approximation of the derivatives 
at vertical walls by one-sided differences. At the horizontal walls, we always use one- 
sided differences. For illustration purposes, in Fig. || we present part of the cosine curve 
evolving under ( |2.5| ) with the speed function F = —0.1 K TT . Boundary conditions 
at vertical walls are periodic. Note that the speed of evolution slows as the curve 
approaches equilibrium state with K = (line y = 0.5). This is because the curvature, 
and hence its derivative, become smaller. In order to demonstrate the abilities of the 
method, in Fig. ||| we present the evolution of a non-smooth curve (step function) under 
the same speed law. 

Remark 4: The very special feature of the presented implementation of the Level 
Set Method is the incorporation of physical boundary conditions into the Level Set 
numerical scheme. Most of the implementations known so far lack this complication. 
Usually only closed interfaces far away from any boundaries domains are considered 
while the evolution proceeds far away from the boundaries. 
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Figure 4: The evolution of a non-smooth curve (step function). The grid used is 100 x 
20000 time steps were made. 



For the GB grooving by surface diffusion, two boundary conditions at the groove 
root are essential: these are conditions of type (2.4), reflecting the fixed slope of the 
interface and the absence of GB atomic flux. The boundary conditions we impose at 
x = l\ are zero slope of the interface and zero flux. The first condition echoes the 
initial flat interface. The second condition guarantees the conservation of matter, i.e. 
a constant area under the groove profile during the evolution. 

Special attention was given to the treatment of these boundary conditions within 
the framework of the Level Set method. Two methods were developed. 

The simplest technique is the use of correction step in the iterative algorithm. 
The fixed slope at the groove root is achieved in the following way: at every time 
step, the interface is reconstructed from the (/>-field and the locations of the two 
end-points of the interface (at x = and x = l\, respectively) are corrected in 
order to preserve the small-slope and the zero-slope conditions. Then, for all grid 
points that lie on grid lines x = and x = h, it is sufficient to directly compute 
a new signed distances to the updated locations of the interface end points. This 
way we incorporate the new locations of the end points back into the 0-field. This 
direct reinitialization is performed only for a few grid points that lie on vertical 
boundaries and, besides, this computation does not contain iteration loop. The 
zero flux conditions could be imposed locally, i.e. in the vicinity of the groove root 
and of the interface end-point at x = l\, or along the the entire x = and x = l\ 
grid lines. After the computed values of K T are reset to zero, the K TT is computed 
according to eq. (|2.10| ), where K T = at x = 0, l\ and K T ^ otherwise. After 
multiplication by — B, this gives the values of the normal velocity function Q2.S| ), 
corrected by the zero flux constraint. 
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Extension of the 0-neld beyond the GB makes use of Taylor expansion up to 
second order, as follows (also see eq. Q2.6| )): 



= 0Oj - <Px \o,j Ax = (p j - \ V(p ,j\ n x |oj Ax = 4> j + |V0 o ,il sin 6» Ax, 

(2.17) 

where 4>—\a is one grid point beyond the GB. Equation ( |2.17j ) incorporates the 
groove root angle. Then we compute in ([2.7D the curvature values, Koj, along the 
GB, using both the values of <j> inside the computational domain (4>ij) and outside 
(4>-i.j)- This also gives us the values of K y |oj . The zero flux condition is applied 
using equation (|l.3| ) which, after substitution of normal vector components from 
( |2.6| ) and rearrangement of the terms become 

K x \ j — \ 0: j — —K y | ,j tan^o- (,^- i °J 

<Py 

Applying Taylor expansion again, we get the ghost values of the curvature: 

K-ij_ = K 0tj - K x |oj Ax, (2.19) 



where K x \qa is given by ( 2.18 ). Now all the data is known and we can compute 
the values of K TT from Q2.10 ) and the values of the normal velocity from ( |2.8[ ). 

Both methods were successfully used in calculations. 



3 Numerical results and discussion 

Figures ^ to ^ show the groove profile having different slopes at the groove root, evolving 
under ( |2.5| ) with a speed function F = —BK TT . We take B = 0.025. The profile is 
symmetric with respect to the GB at x = 0, therefore only its right part is calculated. 
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Figure 5: GB grooving by surface diffusion. The slope at groove root is m = 6.55e — 02. The 
initial interface is shown with dashed-dotted line, the numerical results obtained by means 



of the LS Method are shown with solid lines, the reference results of [14]] are shown with 
dashed lines. 
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Figure 6: Same as Figure BL but the slope at groove root is m = 9.85e — 02. 
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Figure 7: Same as Figures |, §, but the slope at groove root is m = 1.32e — 01. 



The results obtained by means of the LS Method are shown with solid lines, while 
reference results for Mullins' problem (|2.2|)-(2.4) are shown with dashed lines. In all 



the three numerical experiments reported here the dimensions of the computational 
box are [0., 0.08; 0., 0.02], the mesh is 120x40. 

Our initial interface for the Level Set simulations already has the shape of Mullins' 
groove. The reason we don't have a flat interface y(x, 0) = const, as an initial condition 
is that the LS formulation requires a non-zero initial curvature, otherwise the curve 
does not evolve at all (since F = in this case). The initial interface in Figs. [B| - [7| is 
shown with dashed-dotted line. 

The initial Mullins' groove is obtained as follows: we integrate numerically the 
equation (^^) using the method-of-lines approach. The time integrator is second- 
order Runge-Kutta and the spatial operator is discretized using second order central 
differences. The integration proceeds from t = to t = 8.0e — 09. The initial and 
boundary conditions are (|2.3| ) and Q2.4|) , where Oq = 7r/48, 7r/32, 7r/24 stands for Figs. 
|5|- fj], respectively. The corresponding slopes are m = 6.55e — 02, 9.85e — 02, 1.32e — 01. 
The practical values used in experiments lie between 0.05 to 0.2 and the range of the 
groove depth in experiments is between 0.1// and In. The reason we anticipate the 
use of the analytic solution to the Mullins' problem fl2.2|) -(^~4|) (either it exists) is the 
truncation of infinite series in which this solution is represented. The reference results 
for later times are also obtained using the described numerical procedure. 



In 1 14], two kinetic laws were established (within the framework of the SSA). One 
concerns the evolution of the depth of the groove with respect to the maximum surface 
elevation (see Fig. 1). The depth, d, is governed by 



d = 0.973 m (Bt) 



1/4 



(3.1) 



14 



The other kinetic law concerns the evolution of the distance between the position of 
the groove root and that of the surface maximum. In the case of the symmetric groove, 
we call it the half-width, w, of the groove. It is governed by 

w = 2.3(5t) 1/4 . (3.2) 

From these expressions, we have the time independent ratio 

w/d = 2.3515/m. (3.3) 

Under typical experimental conditions a groove of depth d = 0.3/x is formed within 
t = 10 4 sec (2.4 hr). It is shown in |14[ |, that it would require approximately 8 days 
to triple this depth. This explains why in our numerical experiments the groove seems 
to stop developing at later times. The physical reason for this is the increase in the 
length of a path along which the surface diffusion takes place. As a rule, we stop the 
run when the groove doubles its depth or width. 

For the slopes considered, we observe good qualitative agreement with Mullins' 
solution. The small difference is due to two reasons. First, the results to which we 
compare are obtained by integrating the linearized equation ( |2.2| ), which is, strictly 
speaking, valid only for infinitesimal slopes. The slopes we choose are, of course, finite, 



and the governing equation we solve, i.e., the equation (2.5) is fully nonlinear. Second, 
there are inevitable area losses, since the LS method is not fully conservative. For 
bigger slopes, our grooves appear to be deeper and wider than Mullins' one. 
In Tables 1 to 3, the results for all three tests are summarized. 



Table 1. Our results for GB grooving, compared with classical Mullins' results. The 
slope at groove root is m = 6.55e — 02. 



step 


t 


d,eq.(|En) 


d,LS M. 


w,eq-<B 


w,LS M. 


u?/d,eq.(pOD 


w/d,LS M. 





8.0e-9 


2.39e-4 


2.39e-4 


8.60e-3 


8.60e-3 


3.60e+l 


3.60e+l 


2e+3 


1.6e-8 


2.85e-4 


2.50e-4 


1.03e-2 


1.01e-2 


3.60e+l 


4.03e+l 


4e+3 


2.4e-8 


3.15e-4 


2.68e-4 


1.14e-2 


1.08e-2 


3.60e+l 


4.02e+l 


6e+3 


3.2e-8 


3.39e-4 


2.84e-4 


1.22e-2 


1.13e-2 


3.60e+l 


3.99e+l 


8e+3 


4.0e-8 


3.58e-4 


2.99e-4 


1.29e-2 


1.19e-2 


3.60e+l 


3.96e+l 


10e+3 


4.8e-8 


3.75e-4 


3.13e-4 


1.35e-2 


1.23e-2 


3.60e+l 


3.94e+l 


12e+3 


5.6e-8 


3.90e-4 


3.26e-4 


1.41e-2 


1.28e-2 


3.60e+l 


3.91e+l 


14e+3 


6.4e-8 


4.03e-4 


3.38e-4 


1.45e-2 


1.32e-2 


3.60e+l 


3.89e+l 


16e+3 


7.2e-8 


4.15e-4 


3.50e-4 


1.50e-2 


1.35e-2 


3.60e+l 


3.87e+l 


18e+3 


8.0e-8 


4.26e-4 


3.61e-4 


1.54e-2 


1.39e-2 


3.60e+l 


3.85e+l 
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Table 2. Same as Table 1, but the slope at groove root is m = 9.85e — 02. 



step 


t 


d,eq.(|3_l) 


d,LS M. 


w,eq.(3j) 


w,LS M. 


io/d,eq.(3.3|) 


w/d,LS M. 





8.0e-9 


3.59e-4 


3.59e-4 


8.61e-3 


8.61e-3 


2.40e+l 


2.40e+l 


2e+3 


1.6e-8 


4.29e-4 


3.95e-4 


1.03e-2 


1.03e-2 


2.40e+l 


2.61e+l 


4e+3 


2.4e-8 


4.74e-4 


4.38e-4 


1.14e-2 


1.13e-2 


2.40e+l 


2.59e+l 


6e+3 


3.2e-8 


5.10e-4 


4.77e-4 


1.22e-2 


1.21e-2 


2.40e+l 


2.55e+l 


8e+3 


4.0e-8 


5.39e-4 


5.12e-4 


1.30e-2 


1.29e-2 


2.40e+l 


2.52e+l 


10e+3 


4.8e-8 


5.64e-4 


5.45e-4 


1.35e-2 


1.36e-2 


2.40e+l 


2.49e+l 


12e+3 


5.6e-8 


5.86e-4 


5.76e-4 


1.41e-2 


1.42e-2 


2.40e+l 


2.47e+l 


14e+3 


6.4e-8 


6.06e-4 


6.05e-4 


1.45e-2 


1.48e-2 


2.40e+l 


2.44e+l 


16e+3 


7.2e-8 


6.24e-4 


6.33e-4 


1.50e-2 


1.53e-2 


2.40e+l 


2.42e+l 


18e+3 


8.0e-8 


6.41e-4 


6.59e-4 


1.54e-2 


1.58e-2 


2.40e+l 


2.41e+l 



Table 3. Same as Tables 1 and 2, but the slope at groove root is m = 1.32e — 01. 



step 


t 


d,eq.$3) 


d,LS M. 


«j,eq.(|3) 


w,LS M. 


u>/d,eq.(|3.3|) 


w/d,LS M. 





8.0e-9 


4.80e-4 


4.80e-4 


8.61e-3 


8.61e-3 


1.79e+l 


1.79e+l 


2e+3 


1.6e-8 


5.74e-4 


5.60e-4 


1.03e-2 


1.06e-2 


1.79e+l 


1.89e+l 


4e+3 


2.4e-8 


6.36e-4 


6.42e-4 


1.14e-2 


1.19e-2 


1.79e+l 


1.85e+l 


6e+3 


3.2e-8 


6.83e-4 


7.15e-4 


1.22e-2 


1.30e-2 


1.79e+l 


1.81e+l 


8e+3 


4.0e-8 


7.22e-4 


7.80e-4 


1.29e-2 


1.39e-2 


1.79e+l 


1.78e+l 


10e+3 


4.8e-8 


7.56e-4 


8.39e-4 


1.35e-2 


1.47e-2 


1.79e+l 


1.76e+l 


12e+3 


5.6e-8 


7.86e-4 


8.94e-4 


1.41e-2 


1.55e-2 


1.79e+l 


1.74e+l 


14e+3 


6.4e-8 


8.12e-4 


9.44e-4 


1.45e-2 


1.62e-2 


1.79e+l 


1.72e+l 


16e+3 


7.2e-8 


8.36e-4 


9.90e-4 


1.50e-2 


1.69e-2 


1.79e+l 


1.70e+l 


18e+3 


8.0e-8 


8.59e-4 


1.03e-3 


1.54e-2 


1.75e-2 


1.79e+l 


1.69e+l 



An interesting simple extension of the classical two-grain model is the case of a 
periodic array of grains separated by parallel GBs. In Fig. ||, we present the results 
for the evolution of a surface profile intersected by two GBs, i and i + 1. The physical 
boundary conditions at both groove roots are a constant slope of the surface and zero 
flux (for this example, the slope at groove roots is m = 9.85e — 02). At short times, 
grooves develop at each grain boundary according to the solution for an isolated grain 
boundary, as presented in Figs. [5] - [7[ grooving stops when, at sufficiently long times, 
identical circular arcs develop connecting adjacent GBs. The same result was obtained 
in H using Fourier method and the SSA. 
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Figure 8: Long-time evolution of surface profile intersected by two adjacent GBs. The initial 
surface for LS simulations is shown with dashed-dotted line. 



4 Conclusions 

The Level Set method was used to model the grain-boundary grooving by surface 
diffusion in an idealized polygranular interconnect which consists of grains separated 
by parallel GBs. The novel feature of the method is the treatment of physical boundary 
conditions at the groove root. The results obtained are in good agreement with the 
classical one (Mullins, |14|l ) for the case of an isolated grain boundary (two-grain case) 
and with more recent results of Q for the case of periodic array of grains. One goal for 
future work is to apply electromigration influence on the grooving process. In addition, 
the algorithm and its software implementation will be used by materials scientists to 
pursue studies of GB grooving with an arbitrary electromigration flux, the various 
ratio of the GB to surface diffusivity which was predicted to critically affect the groove 
kinetics and shape account for various EM failure regimes 
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